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We present predictions for the temperature dependent shifts and damping rates. They are obtained 
by applying the dielectric formalism to a simple model of a trapped Bose gas. Within the framework 
of the model we use lowest order perturbation theory to determine the first order correction to the 
results of Hartree-Fock-Bogoliubov-Popov theory for the complex collective excitation frequencies, 
and present numerical results for the temperature dependence of the damping rates and the frequency 
shifts. Good agreement with the experimental values measured at JILA are found for the m = 2 
mode, while we find disagreements in the shifts for m — 0. The latter point to the necessity of a 
non-perturbative treatment for an explanation of the temperature-dependence of the m=0 shifts. 
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I. INTRODUCTION 

Since the discovery of Bose-Einstein condensation in 
traps a wealth of experimental data on coUective excita- 
tions has appeared in the hterature (for experimental re- 
views see 1^,^) waiting for theoretical explanation. Some 
of the earliest measurements, still requiring a firm theo- 
retical analysis, were performed in oscillating traps which 
permitted the selective excitation of different excitation- 
modes 1^-^ . Common features of the above experiments 
are that the density of atoms in the trap is relatively 
small and the temperature is extremely low. Conse- 
quently, from the theoretical point of view the atoms 
can be treated as a weakly interacting degenerate Bose 
gas, while the interaction potential is well described by 
the s-wave approximation. Due to the presence of the 
trap the whole system is not translationally invariant, 
and field equations of any approximation must be solved 
in real-space, not in momentum space. Furthermore the 
excitation spectrum is discrete rather than continuous as 
in spatially homogeneous condensates like helium II. 

Nevertheless, most theoretical approaches are based on 
the natural generalization of one or the other of the ho- 
mogeneous descriptions to the inhomogeneous case. The 
present paper also belongs to this line. It is based on the 
dielectric formalism (see Ref. Q and further references 
therein), first introduced for spatially homogeneous sys- 
tems at zero temperature [p|-p^ , later used at finite tem- 
perature and recently generalized to inhomoge- 
neous systems in |13 . The great success of the dielectric 
formalism lies in showing that the order parameter cor- 
relation function (the one-particle Green's function) and 
the density-density correlation function have the same 
spectra below the critical temperature. In principle the 
dielectric formalism is valid at all temperatures but still 



requires to deal with infinitely many graphs to obtain 
exact results. In practice however, one resorts to some 
approximations for the proper and irreducible part of cer- 
tain quantities of physical interest, which are the key 
quantities in the dielectric formalism. Following and ex- 
tending Ref. ||ll|,|l3| we shall base our investigation on a 
simple model for trapped weakly interacting Bose gases, 
which in homogeneous systems (a) is valid at finite tem- 
perature, (b) satisfies the generalized Ward identities, 
and (c) guarantees that both the Green's function and 
the density-density autocorrelation spectra exhibit the 
same excitations. In addition the model we use, which 
builds on and extends a simpler one discussed in detail 
in Ref. , can be shown to be related to an approxi- 
mation used by Minguzzi and Tosi [ p^ . 

It is our main purpose here to evaluate, within the 
approximation defined by our choice of the model, the 
damping rates and frequency shifts and to compare the 
results with the experimental data Q of the JILA group. 
In principle we are also able to calculate theoretical val- 
ues for the MIT measurement, but in our present ap- 
proach we treat the fluctuations of the thermal den- 
sity perturbatively and consider only Landau damping. 
This approach is justified if the condition fc^T 3> is 
fulfilled which is only the case for the JILA measure- 
ments. Furthermore, the high anisotropy of the MIT 
trap leads to some numerical difficulties in the code we 
use so far. ^From the theoretical side several papers 
have appeared in the literature going beyond the Hartree- 
Fock-Bogoliubov-Popov (HFBP) theory which is re- 
ally necessary to take into account damping processes. 
Most relevant papers beyond HFBP describe the damp- 
ing process p^-|20[| or calculate the shifts by including the 
anomalous average in such a way that the resulting ap- 
proximation is gapless [ pl| . The approach applied in Ref. 
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pO| is the second order Behaev theory which is known to 
be gapless. It treats both Behaev and Landau dampings 
and also calculates the shift of elementary excitations in 
local density approximation along with their damping. 

In the present paper we wish to present a theory of 
the shifts and widths of the low-lying modes. Our ap- 
proach is based on the dielectric formalism and achieves 
its simplicity by a judicious selection of a subset of graphs 
for the proper part of the physical quantities of interest. 
The model we shall investigate in detail in this work ac- 
counts for the Landau damping, but cannot account for 
Beliaev damping. This sacrifice for gain in simplicity is 
not too big because in the temperature region where the 
measurements we wish to understand are performed, and 
where we calculate the shifts and the damping rates, the 
Beliaev damping of the excitations is negligible. 

The paper is organized as follows. In Sec. || we 
briefly summarize the general framework of the dielec- 
tric formalism for inhomogeneous systems and consider 
an extended version of the model approximation of Ref. 
1^,^. The conditions for choosing the necessary ba- 
sic building blocks for the ladder approximation and for 
the proper and irreducible quantities are spelled out. 
Then the quantities given in our formalism are related to 
the fluctuations of the condensate and the thermal den- 
sity. We show that neglecting the thermal density fluc- 
tuations we recover the usual Hartree-Fock-Bogoliubov- 
Popov equations ||l^. Then we derive the corrections of 
the HFBP excitation energies to first order in the thermal 
density fluctuations. To solve the closed set of equations 
for the damping and the shifts still requires some numer- 
ical work in the inhomogeneous case. We briefly discuss 
our numerical procedure in Sec. 



III. Sec. [V 



contains 

the discussion of our results and the comparison with the 
experimental data measured at JILA. The last Section 
is devoted to the conclusions and to some final remarks. 
There we also compare our model approximation with 
other approximations given in the literature, e.g. the al- 
ready mentioned treatment of all the Beliaev diagrams 
by Fedichev et al. , the kinetic equations of Ref. |lj] 
and the coUisionless Bolzmann equation |23,M. 



II. FORMULATION 

Here we summarize the dielectric formalism first ap- 
plied to inhomogeneous systems in [ p^ . However, we 
shall not repeat the whole treatment, but rather con- 
centrate on the key points, and indicate the new further 
steps. 

The Hamiltonian of our problem in second quantized 
form is 

H = J d^r *t(r) (^-^ A + C/(r)^ «'(r) 

+ i / (fri / d^r2 *^(ri)*^(r2)-y(ri,r2)#(ri)*(r2), 



(1) 

where ^'(r) is the Bose field operator, U{r) is the trap 
potential, and v{ri,r2) describes the two-body interac- 
tion. In the following it is chosen as 



v{r, r') — gS{r — r') 



S{r ~ r') 



(2) 



where a is the s-wave scattering length and m is the mass 
of the atoms. Throughout we shall restrict ourselves to 
temperatures below the critical temperature. As usual, 
for T < Tc we split off the condensate wave-function 
<fo(r), 

*(r) = $oW + $(r), (3) 

where <I'o(''') = {''')) and (. . .) denotes thermal averag- 
ing 



Trie-'3(g-M^) 
Xre-/9(^-M^) 



(4) 



Since we are interested in finite temperature excitations 
we use Green's functions 



Ga,0{r,T;r',T') 



(5) 



with field operators in Matsubara representation 
$i(r,r) = $(r,r) and $2(^",t) = ^^{r.r). The other 
key quantity we are interested in is the density autocor- 
relation function 



X(r, r; r', r') - -(T, [n(r, r)n(r', r')]), 



(6) 



where fi(r) — n{r) — {n{r)). The finite temperature 
Green's functions (||) and the autocorrelation function 
are functions of r, t' via r — r' only and are periodic 
with period (3h. Thus, one can expand them as Fourier 
series. The Matsubara Fourier coefficients are given by 



/3h 



ciTe'-"^G„,^(r,r',r), 





2n7r 



(7) 



where n is an integer. A corresponding expansion is made 
for x(r, r; r', r') with coefficients x(r, r', iuin)- Retarded 
functions can be obtained by the usual analytic continu- 
ation (iujn —^uj + irj, rj infinitesimal). 

It is useful to introduce the proper part of a quantity 
which is defined as the sum of diagrammatic contribu- 
tions, which cannot he split into two parts by cutting a 
single interaction line. In the following proper parts will 
be denoted by a tilde. By definition the density autocor- 
relation function and its proper part fulfill 



X(r,r',w) = x{r,r\uj) 

xx(r2,r',tj). 



d^r2x{r,ri,uj)v{ri,r2) 



(8) 
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When x(r,r',a;) has a pole in but its proper part is 
nonsingular at the same w, then there exists an eigen- 
function ^(r) satisfying 

?(r) = ^J dVi J d^r2xir,ri,uj)v{ri,r2Mr2). (9) 

This eigenfunction can be identified with the total den- 
sity fluctuation Sn{r) = ^{r) at the eigenfrequency given 
by the pole. In a similar way we get the eigenfunctions 
(/?i,(^2 corresponding to the 1-particle Green's functions 
Gq,/3 by solving the eigenproblem: 

<Pa(r) = ^J (fri J (fr2G^^{r,ri,uj)^^^f}{ri,r2,uj) 

X'^/3(r2) (10) 

(Here and in the following we adopt the convention that 
summation has to be taken over repeated greek indices). 

GfO(r,r',c^„) = GrO(r',r,-c.„) 

^f"(r)^fQ(rO* 



Gg^{r,r',u;r,) = Gg^ir,r',u;,,)^0 

are the Green's functions of the free harmonic trap with 
the corresponding harmonic oscillator excitation frequen- 
cies ef^, eigenfunctions ip^'^{r) and chemical potential 
jjL^'-' = {Tl/2){u!j. + u!y + oJz), where ojx.y,z Siie the three 
trap frequencies. Sq,/3 are the selfenergies describing the 
corrections due to the interactions. The condensate den- 
sity fluctuations Sndr) given by perturbations of the con- 
densate density ndr) = |$o('")P due to excitations of 
single-quasiparticle modes can be identified with 



(11) 



Below Tc the proper part x can be further decom- 
posed into irreducibile (also termed "regular") x^^'^ and 
reducible (also termed "singular") x^*^ parts. We call a 
diagram irreducible if it cannot be split into two parts by 
cutting a single particle line. The reducible part x^'^'^ is 
related to the existence of the so-called anomalous proper 
vertex Aq, which contains all the proper diagrams with 
only one outer interaction line and only one outer particle 
line, and which is due to the presence of the condensate 
below Tc'. 

X{r,r',u;) = x^'''> {r , r' , uj) + x^''> {r , r' , uj) 
X^'')(r,r',u;) = J d^ri j dVa Aa(r, ri, w)G„,/5(r-i, ra, w) 
xA^(r2,r',Lj) . 

(12) 

Here Ga,i3{i'i,r2Tijo) is the proper part of the Green's 
function satisfying Dyson's equation with the proper part 
5^^q,/3(t"i, r2, w) of the selfenergy only. In Ref. ||l^] it 



is shown that the eigenfunctions ^, ipa belonging to the 
same eigenvalue w are related by the anomalous vertex 
Aq containing all the proper and improper vertex contri- 
butions 



i{r) = J d^riAa{r,ri,uj)(pairi) 



(13) 



So far no approximation has been made. Now we define 
approximate expressions for the building blocks A'^ and 
X° of the ladder approximation in which we derive the 
proper quantities A^ and x*^*^'- 
The trivial vertex function 

A;(r, r', L,) = AO(r, r' , u) = $o(r)(5(r - r') (14) 

and the bubble graph 

X°(r,rV) = 



= - ;4 E f (r', r, zc.„ - c.) (15) 



with the Green's function G^^ corresponding to the 
Hartree-Fock Hamiltonian 



i7«^(r) 



2V72 



2m 



+ Uir) + 2gn{r). (16) 



The Green's bmctions G^^ satisfy the relations: 

{huj-H"^{r))G^/{r,r',uj) = hS{r-r') , (17) 
G2Y(r,r',c^) = Gf/(r',r,-c.) , (18) 
Gf2^ = G2T = 0. (19) 

We have already used the Hartree-Fock Green's functions 
in the bubble graph. We will see later that this approach 
turns out to be consistent. 

In the ladder approximation the proper contributions 
X^'^\r , r' , Lo) are derived by subsequent insertions of in- 
teraction lines into the bubble diagrams. Therefore 
^/^ jg determined by the selfconsistent equation: 

x'-''\r,r',uj) = x°ir,r',uj) 

+ 1 J dVx°(r,ri,c.)x('-)(ri,r',c.) (20) 

For the proper vertex function A^ we take into account 
the trivial vertex A[^(r) and the first order correction 

Aa\r,r',u!) = gx^ifTr' ,u!)^o{'''')- In A^^ we replace 
the single interaction line by the T-matrix which defines 
the following equation for the proper vertex function Aq : 

A4r,r',c.)=A«(r,r',c.) 

+ 1 / dViX°(r,ri,c^)A„(ri,r',c^), (21) 
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from which it can be determined selfconsistently. 

The analogous Dyson equation for the complete vertex 

function given by h.a{r,r',bj) — 

Aa{r,r',uj) + f / d^ri x'''^\r,ri,uj)Aa{ri,r',uj) can be 
easily expressed in terms of the building blocks of the 
ladder approximation: 

A„(r,r',c.)=AO(r,r',c.) 

+2| I dVx°(r,ri,c.)A„(ri,r',c.) (22) 

We identify the thermal density fluctuations Snrir) = 
dn{r) — Snc{r) by inserting ( ]2^ ) and ( pi] ) in (13) and 
using (|ll|) and (|l|): 

SnTir,uj) ^2^ J d^rix°{r,ri,uj)Sn{ri,uj) (23) 

Another way to derive this result is to consider the pos- 
sible diagrams in the linear response function xt for the 
thermal density where xt is defined by: 



Snrir)^ (fr'xT{r,r',uj)SV{r') 



(24) 



and SV is an additional small perturbation coupling to 
the density operator. 

The only diagrams of our approximation for x not con- 
tributing to XT are those which start with the trivial ver- 
tex function A° — $o on the side coupling to the thermal 
density fluctuation: 

XT{r,r\uj) = j (fri xir,ri,uj) - J d^r2 J d^rs 

A" (r, r2, uj)Ga^f3{r2,r3, w)A^(r3, ri, w) 
X 5{ri-r') + ^xirur',Lu)\ (25) 
= J d^ri x^''\r,ri,uj) 

+ 1 J d^r2x"{r,r2,uj)x^''>{r2,ri,uj) 
5{ri-r') + ^x{rur',iu)\ (26) 

In the previous step we have inserted the expression for 
A° given by equation (|2|) into (|25|). If we use 6V = 
J x"^'^"' we get 

Snxir) ~ / d'^ri / d'^r'xT{''',ri,LL!)x~^{ri,r' ,uj)5n{r') 



d^ri I d^r2 I d^r' 



X^'^\r,ri,uj) 



+ -X°ir,r2,uj)x^'Hr2,ri,uj) 

X-Hri,r', uj) + f 5(ri - r')] 6n{r') (27) 

which reduces in the case of resonance (J x^^'^"- — 0) to 
the relation: 



Snrir) = / d^ri / d^r' \x^''\r,ri,u) 



Using Sn - 



+ ^x'{r,r',u;)x(^Hr,r',Lu)] |<5n(r') (28) 
Sric + SriT the equation for duc has the form: 
= J d'r, J d^r'{6{r-r,)-^x°{r,r,,uj)) 



xx^''^(ri,r',w)|(5n(r') 



(29) 



^From equations (28) with (j2^) and (||) it is straightfor- 
ward to derive the final result Snx = 2{g/h) J x^Sn. 

For the purpose of calculating Sric we need to make ad- 
ditional approximations for the proper selfenergies Sa./j. 
We restrict ourselves to the proper and irreducible di- 
agrams which (a) are only 0-loop and 1-loop diagrams 
according to the classification by Wong and Gould |l^] 
(justified by assuming Sut to be small and only calculat- 
ing its first order contributions), (b) are connected with 
Landau damping (since wc neglect the Beliaev damp- 
ing) and (c) do not contain anomalous Greensfunctions 
(Popov approximation). 

First we introduce those selfenergies ^ appearing in 
the Gross-Pitaevskii-equation for the condensate: 



HoMr) = 0. 
where Hq is given by 



(30) 



Ho = - — A + U{r) - /i + 5l*o(r)P + 2.gnT(r), (31) 
2m 

t%{ryM = {9mr)\'' + 2gnr{r,))[l \^ (32) 

They contain the 0-loop diagram g|<I'o('")P together with 
the contributions from the stationary density of the non- 
condensed atoms: 



nT{r) = (l>'''(r)l>(r)). 



(33) 



This density can be calculated in two different ways: ei- 
ther in Popov approximation as in [p^ , which is a self- 
consistent, gapless approximation, but does not satisfy 
the Ward identities in the homogeneous case, or by us- 
ing the Hartree-Fock approximation, which neglects the 
quasi-particle aspects of the thermal excitations but has 
the great virtue of satisfying the Ward identities in the 
spatially homogeneous case Therefore we choose 

this second method in the following. However, it is re- 
markable that the final results for the shifts and widths 
we obtain are very similar for both methods of calcula- 
tion. Due to conditions (a)-(c) we are left with only four 
additional selfenergy contributions which, using the same 
T-matrix approximation as before, can be written as: 



sW(r,r',c^) = Mr-,r-',^)( { { 



1 1 



(34) 
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with 



The corresponding Greensfunctions 

^1.2 — ^2 1 — 0, 



(35) 



(36) 



Gl,{r,r',u;)^Gl,{r',r,-co) (37) 

satisfy 

{hLu ~ i/o)G° i(r, r', to) = M(r - r'), (38) 
and ^a/3(''7 ''''i ^) is determined by: 
J2Gc.p{r,r',u;)^ J d'r, J d^r^ ^ GO,(r, n, c.) 

x((5(ri - r2)(5(r2 - r') + a{ri,r2,uj) 
xY,Gap{r2,r',Lo)) . (39) 



Q/3 



Inserting x'*^ = X]q/3 ^aGa/jA^ in equation (29) we get 
for dnc{r): 

9 10, 



Ax{ri,r2,Lu)Sn{r2) + |A;((ri, rg, w) 



(40) 



xx^'''>{r3,r2,uj)Snc{r2) 
which may be rewritten as 

Sn,{r) = J d^r, J d3r2|<i>o(r) ^ G°^(r, n, c.)<i>o(ri) 



Q/3 



5{ri - r2)Sn{r2) + ^X^''\ri,r2,uj) 
X ((5n(r2) + 6nc{r2))] . 



(41) 



To solve Eq. (|41j) is still ver y d ifficult. Instead, we fol- 
low the procedure applied in [|l3[, where the term f x*''^'' 
was treated as a small perturbation and the excitation 
frequencies were determined using first order perturba- 
tion theory. In the perturbation term j-X^^^ we may re- 
place x^^^ by and Suc by Sn to first order. We then 
use eq.(23) to obtain 



Sn,{r)^ [ d3ri|$o(r)^G°^(r,ri,c.)<i>o(ri) 

X [5n,{r^) + 26nT{r^)] . (42) 
In the next step we show that the unperturbed problem 

<5n°(r) = | / rf3ri<i>o(r)^G°^(r,ri,c.)$o(ri) 

•' a/3 



X(5n°(ri) 



(43) 



is equivalent to the Hartree-Fock-Bogoliubov-Popov cal- 
culation of 5n^. By dividing equation (|4^ ) on both sides 
with the condensate function <I>o('') and afterwards muli- 
plying both sides from the left with {huj — Ho){—fiuj — Ho) 
we derive an equation of the form of the diagonalized 
HFBP equations: 



$o(r) 



-2gHo{r)nc{r) 



$o(t-) *o(r-) 



(44) 



This equivalence can be shown by performing the sum 
and the difference of the Hartree-Fock-Bogoliubov equa- 
tions: 



gnc{r) H°{r)+gnc{r) 



(45) 



We get the following diagonalized equations for lyj^-''' ± 



h'4'\^^^\r) + ^'i\r)) = [Hl{r) + 2gffo(r)n.(r) 

x{v^f{r) + v^i\r)) (46) 
Hl{r) + 2gn,[r)Ho{r) 



So), 



xi^['\r)~ip'^'{r)). (47) 



0)/ 



If we set SriT — (unperturbed case) in equation ( |44| ) 
we recover the HFBP equation where corresponds to 

the eigenvalues — and Sn^ to the corresponding 

eigenf unctions Sn° = — ^oiVi^ + '^2"'^)- 
The first order correction due to the dynamics of the ther- 
mal density Snx (or in other words due to the existence 
of x") can be calculated by expanding 



Sric 



■ 1^1 + 



(48) 
(49) 



where 6nl and uji are the first order corrections. 
Inserting this perturbation ansatz in ( ^ ) the first order 
contributions must satisfy the equation: 



2, .2 



H^ir) + 2gHoir)n,{r)-h 



2nVc^o$^-45i?o(r)n,(r) 



$o(r) 



Mr) 
SnT{r) 
Mr) 



(50) 



Multiplying this equation from the left with {ipi{r) — 
'P2{r)) and integrating over space the left side of the 
equation vanishes. Due to the normalization condition 
J d^r{ipi{r) — ip2{r))((pi{r) + (p2{r)) = 1 and the rela- 
tion Ho{ipi — <y52) = ^'^o(vi V2) we get the following 
result for lui: 
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huJi = 2g J (frSnl{r)SnT{r) 



(51) 



= d^'r' dnl{r)x°{r,r' ,uo)Sn,{r'). (52) 

This is one of the basic resuhs used in the following for 
the calculation of shifts and widths of collective excita- 
tions. The consistency condition for the applicability of 
perturbation theory is that loq ^ \uJi\ must hold, which 
we shall check below. 

Next, let us concentrate on the quantity {r , r' , ujo) . 
If one expresses G^f (r, r', iujn) in terms of the eigenval- 
ues and eigenvectors of the Hartree-Fock operator H^^: 



then the standard Matsubara sum in Eq. ( p^ ) can be 
easily performed 



X°(r,r',zc.„) = V(n(£f^)-n(e 



vr {r)^r*{r')^f'ir')vr*{r) 



HF 



-HF 



(54) 



where n{e^ ) is the Bose-factor. ^From the form of the 
denominator it is clear that we treat processes leading to 
Landau damping, but not to Beliaev damping. In calcu- 
lating the retarded function via the analytic continua- 
tion iiOn uo + ir] its imaginary part contains a series of 
Dirac-delta peaks due to the discrete spectrum of H^^ . 
In order to get finite damping (i.e., finite imaginary part 
of wi in Eq. (|5^)), which was measured in experiments, 
one needs some smoothing over the Dirac-delta peaks. 
In principle, for the shift there is no need for smoothing, 
but one must use the same procedure both for the real 
and the imaginary part of , otherwise they do not ful- 
fill the Kramers-Kronig relations. The smoothing is done 
by using the semiclassical or local density approximation. 
First, we calculate the semiclassical Green's function. To 
do this, let us consider the representation of the Green's 
function in the center-of-mass, k representation, defined 
as: 



Gfr(r,r',. 



(Pk 
(2^ 



-r')l2,k,u) 
(55) 



For this representation the semiclassical Green's function 
is simply 



Gff(^^)(H,fc,z^„) 



ihuon - (efc - m(-R)) ' 



(56) 



where 



2m 



H{R) = M - U{R) - 2gnc{R) 



2gnT{R). 
(57) 



It is easy to show that if one uses the mixed represen- 
tation for x^ , similarly defined as in (pq) 



x"{r,r',uj) = 



(Pk 
(2^ 



ik{r- 



'■')x°((r + r')/2,fc,c.) 



(58) 



then in the new representation Eq.([l5|) takes the form 

d^k' 



xGi i {R, k' , iujp — iujn)- 



(59) 



Performing the Matsubara sum, but now with the semi- 
classical Green's function (|56|), after the analytic contin- 
uation we arrive at 



x'{R,k,u;) = 



d^k' 



nk' — rik+k' 



(27r)3 uj - [ek+k' - f-k')lft' + i-q ' 



(60) 



where is the Bose-factor appearing in the semiclassical 
(or local-density) approximation: 



nk 



1 



A' 



(61) 



1 



In the region where the condensate density vanishes the 
condition /i(i2) ~ holds. As a result, in our semiclassi- 
cal approximation we would there get a logarithmic sin- 
gularity of x^{R,k,w) at k = \k\ = \/2mhL0 leading to 
large contributions from frequencies around zero, which 
are completely artificial, since the trapped system has 
no levels below a certain lowest-lying level. We there- 
fore suppress these unphysical logarithmically singular 
contributions by choosing the lowest-lying energy level 
^mm = 2.12a;p of the Hartree-Fock operator as a natural 
lower energy cut-off in our application of the Hartree- 
Fock and local density approximation. Here Wp is the 
smallest trap frequency. We thus ensure that all the con- 
tributions to frequency shifts and damping rates come 
only from interactions with modes in the physical energy 
region above this lowest Hartree-Fock energy level. In 
Fig. 1^ we plot the spectral contributions to the ther- 
mally excited atoms as a function of the frequency for 
both, the local density approximation and the histogram 
of energy-levels obtained by the direct diagonalization of 
the Hartree-Fock Hamiltonian. The agreement above our 
cut-off proves to be rather good so that we can expect 
to get reliable results in the local density approximation. 
With this cut-off the imaginary part of x° becomes: 
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In 



1 — exp 



(ri^+£fc)^ 

4efc 



(62) 



for - 1)2 > max [2mn~'^{n{R) + hLu„,in),0] - In the 
opposite case we have due to the cut-off the fohowing 
expression: 



X hi 



1 — exp {—Phujri 



1 - exp {-(3h {uJrnin + w)) 



(63) 



for - 1)2 < max [2mh^^ {im{R) + ;iw„i„),0]. After 
performing the principal value integral the real part can 
be written as 



k' dk' np. 



Aw^hq Jn^2me(R) 

'a{k,k',uj,R)b{k,k',uj,R) 



X In 



c{k, k' , uj)d{k, k' , lu) 



(64) 



"fc2 


kk' e(R) 




2m 


m ' h? 


2m 


'kk' 




<R)] 


m 


2m ' 2m 





— h ^uj 



with 

a(fc, k' ,uj, R) — max 
6(fc, A:', cij, i?) = min 



c(fc, k — \ n uj 

2m m, 

fc2 kk' 

d{k, k' , u}) = —- h^^uj 

2m m 

e{R) is an abbreviation for 
e{R) = max [^liR) + fiujmin, 0]. 



III. NUMERICAL METHODS 

In our perturbative calculation of loi we should solve 
first the unperturbed problem for ljq and ^^^(t"). This 
step requires to solve the Gross-Pitaevskii (GP) and the 
HFBP equations. We already have reported on our self- 
consistent algorithm in |24| for this problem, which solves 
the GPE for ^^{r) and for the chemical potential us- 
ing multigrid methods, while the HFBP equations are 
solved in a large truncated basis. As a result, we have all 
the necessary inputs, namely <&o('"), + '^2(1'), wo, 

nxir) for performing the second step: the use of Eq. ( |5^ ) 
in connection with the retarded density auto-correlation 
function Eqs. (||) and (|60|). 

This step is straightforward, but from the nu- 
merical point of view it requires to cope with 
six-dimensional space-integrals and three-dimensional 



momentum-integrals. To reduce the number of numer- 
ical integrations we proceed as follows. If we introduce 
C as an intermediate quantity for 



C[R,k) 



d^re-'^''5nl{R- -)5n^{R+ -) (65) 



then Eq. (^2|) can be written as 

?i^i=4^ j d'R J ^x\R,-k,uo)C{R,k). (66) 

Next, let us expand 5nc{r) — $o(^)('i5i(r) + ^2{i')) in 
some isotropic harmonic oscillator basis as 

5n,{r) = <i>o(r)(^i(r) + ^2{r)) = Y,A,U{r), (67) 

where the set of basis functions fi{r) have the same an- 
gular momentum quantum number and parity with re- 
spect to reflection on the x-y plane as the unperturbed 
elementary excitation in question. The next task is to 
find a connection formula, which transforms any prod- 
uct f* {r)fj{r') to the mixed {R,k) representation. This 
step requires some straighforward but rather lengthy cal- 
culation, but the advantage we gain is that due to the 
facts that x^''^(i?, fc, wq) is axially symmetric in R for 
axially symmetric harmonic oscillator trap potential and 
isotropic in k we can perform the integrations over the 
azimuthal angle of R and over the angles of k in (|6^). 
At the end, apart from discrete sums, we have two spa- 
tial integrals in the radial and in the z directions of R 
and one momentum integral for the imaginary part and 
two subsequent momentum integrals for the real part of 
uji. Of course, the remaining integrals must be performed 
numerically. 

The form of the above mentioned connection formula 
which we actually used is based on the following facts: 
First, the product f*{i')fj{r') is a solution of the two- 
body, non-interacting harmonic oscillator problem, which 
is also separable in center-of-mass and relative coordi- 
nates. Thus, expanding it in the other basis containing 
products of harmonic oscillator wave- functions depend- 
ing on R = [r + r')/2 and r — r' the expansion will 
contain only a finite number of new basis functions from 
the energy shell £i + ej . Secondly, the Fourier-transform 
of the harmonic oscillator wave- functions (which we must 
perform with respect to the relative coordinates) are ba- 
sically the same as in real space. 



IV. RESULTS 

We present our numerical results in a way which is di- 
rectly comparable with the measurement [Q. In Figs. ^ 
and we put on the horizontal axes T', which is propor- 
tional to the RF frequency used in evaporative cooling 
and is also the same quantity, used in Fig. 1 of Ref. 
0. From that figure we take the same T = T{T') and 
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A^o = No{T'), which were determined experimentally. 
Also the same anisotropy uJz/ujo as in the experiment was 
used. We present our results as a function of the exper- 
imental T' because different experimental points belong 
to different physical conditions since neither N nor Nq 
could be kept fixed in the experiment. 

In Fig. H we plot the imaginary part of uji together 
with the measured damping rates for the m = 0, 2 modes. 
For temperatures T > 0.6Tc where the Landau damping 
arising from the bubble graph x°(r, r',uj) is the dominant 
process our results are inside the error bars. At lower 
temperatures the calculated damping rates are slightly 
above the measured ones in agreement with the predic- 
tion of calculations for the homogeneous system. There 
the values for the damping rates due to our model calcu- 
lation are larger with a factor of ||- compared with the 
treatment of all Beliaev graphs to second order. This 
is not surprising, because the neglection of the Beliaev 
damping and the use of local density approximation is 
only justified for temperatures satisfying the condition 
ksT » /Lt, only valid for T ^ 0.6Tc. The same condition 
is necessary to neglect the quasi-particle character of the 
additional excited particles contributing to the selfener- 
gies via scattering processes. Comparing our results for 
the damping rates with the values of the unperturbed ujq 
it is manifest that 3wi <ti loq. The real part of uji is 
also much smaller then loq for all the calculated points. 
Therefore the criterium for the applicability of perturba- 
tion theory is justified. 

In Fig. H we plot ujq + SRwi together with the experi- 
mentally measured frequencies. For m = 2 we have good 
agreement for temperatures above T > O.GTc and slightly 
larger negative shifts compared to the experimental val- 
ues below T — .6Tc. The results for the m — mode 
show the same slight difference in the low temperature 
region. But approaching higher temperatures the exper- 
imental and the numerical results predict a different be- 
havior for the excitation energies. While experimentally 
they are increasing we found decreasing values in our 
model calculation. The same qualitative discrepancy has 
also been found by Hutchinson et al. in using a dif- 
ferent approach to the problem. They also reported good 
agreement with the measured frequencies for m = 2, but 
disagreement for m = 0. A possible explanation for the 
above discrepancies has been given by Bijlsma and Stoof 
in Ref. who used the quantum Boltzmann-equation 
in the coUisionless regime and calculated the low-lying ex- 
citation frequencies using a variational approach for the 
TO = 0, 2 modes. They found that there are two nearby 
TO = modes, one of them describing the in-phase the 
other one the out-of-phase oscillations of the condensate 
and the thermal cloud. They concluded that it might be 
possible that in the experiment both are excited together 
and that in Ref. Q only the upper one or a superposi- 
tion of both might be plotted. If there are two nearby 
TO = modes we must explain why in our calculation 
we have obtained only one of them (and furthermore not 
the experimentally measured one). It is clear that by 



our perturbative treatment we fixed ourselves to one of 
them by calculating the corrections to the Hartree-Fock- 
Bogoliubov-Popov equations for Sut ^ Suc and we fol- 
lowed it in changing the temperature. Our results corre- 
spond to the branch describing the out-of-phase motion 
calculated in Ref. [^,^. To see this we recall our for- 
mula for the corrections uji = 2g J <Pr 5n*^{r)5nT{r) and 
the fact that loi is found to be negative in our results. 
Therefore the spacial average over the amplitude prod- 
uct of SriT and Sric is negative which means that Stit and 
5nc are out-of-phase. 

The in-phase mode can accordingly be obtained, if, in 
fact, there is another nearby to = mode in our model 
approximation, if we do not use perturbation theory for 
solving (^^. Instead we should consider the eigenval- 
ues of the bubble graph corresponding to the thermal 
density fluctuation Stit and calculate their perturbation 
due to Suc- This possibility is supported by the early 
study of the behavior of the poles of the density-density 
autocorrelation function in the homogeneous Bose gas 
Iprf , where it was shown that there are more, but purely 
damped modes. They remain purely overdamped above 
Tc in accord with a recent experiment p6t . It can be 
shown that one of the overdamped damped modes cor- 
responds to SriT S> and that it is a mode where Sut 
and Sric relax in phase ||2^ . In addition it is clear that in 
a harmonic trapping potential these damped modes will 
acquire a nonzero real part of their oscillation frequency. 
However to really demonstrate this for the cigensolutions 
of the density-density autocorrelation function in the in- 
homogeneous case, and furthermore to show that in this 
case the missing mode can be described in this way re- 
quires further analytical and numerical studies. 

V. CONCLUSIONS AND FURTHER REMARKS 

In this paper we have investigated the predictions of 
the dielectric formalism applied to the model of a Bose 
gas obtained by extending simpler models discussed in 
the literature |0,^,|3| . The proper and irreducible build- 
ing blocks X*-''-' and A are already the result of summing 
up higher order contributions in the form of a ladder 
approximation only summing up geometric series of the 
bubble graph . ^From the second order Beliaev dia- 
grams we keep for the irreducible selfenergies only those 
which are supposed to give large contributions. Since 
Landau damping dominates Beliaev damping in the tem- 
perature region T ^ 0.6Tc of the measurement at JILA 
we neglect the diagrams connected with Beliaev damp- 
ing. Following the classification of Wong and Gould [lC| ] 
we reduce the amount of diagrams further to only 0- 
loop and 1-loop diagrams. In addition we neglect all the 
selfenergy diagrams containing anomalous Green's func- 
tions. To find poles of the density-density correlation 
function, i.e., collective excitations, we used first order 
perturbation theory, where the proper and regular part 
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of the density-density correlation function played the role 
of the perturbing operator. We showed that the physical 
content of the unperturbed problem is equivalent to find 
elementary excitations in the Hartrce-Fock-Bogoliubov- 
Popov approximation. The fact that all the first order 
corrections to the Hartree-Fock-Bogoliubov-Popov equa- 
tions can be described with one bubble graph permits 
to draw conclusions about the behavior of the quantities 
6nc and SriT from our numerical results. Due to the nega- 
tive energy shifts obtained in our numerical results we can 
identify their motion to be out-of-phase. Since we estab- 
lish our model on the basis of the dielectric formalism we 
achieve the required correspondence of the spectra of the 
Green's functions and the density auto-correlation func- 
tion, ensuring gapless spectra of quasi-particle modes and 
density modes in the homogeneous limit. The approxi- 
mate version (^2|) of our equation (^l|) agrees with the dy- 
namic Hartree-Fock-Bogoliubov-Popov theory presented 
by Minguzzi and Tosi . They derive (p^ ) by lineariz- 
ing the time-dependent Gross-Pitaevskii equation around 
its stationary solution n(r, t) = ^^(r) -|-(5nc(r, t) and ( |2^ ) 
by an equivalent linearization of the Hartree-Fock equa- 
tion for the non-condensate field operator. They show 
that the correct bubble graph x° connected with ( ^ ) is 
a product of Green's functions in Hartree-Fock approxi- 
mation. 

Bijlsma and Stoof |2^] also consider the equation ( p3| ) and 
use the dispersion relation of the Hartree-Fock operator 
in local density approximation for their calculations. In 
contrast to us they examine the coUisionless Bolzmann 
equation for the Wigner function and solve the equations 
with a variational approach, which cannot describe any 
damping of the excitations. Whereas the local density 
approximation we and use is better for large con- 
densates the approach of Fedichev et al. is better 
suited for the opposite case of small condensates due to 
the assumption of hard chaos, which can only be fulfilled 
for small condensates. In this approach the diagrams are 
evaluated by integrating along the classical trajectories. 

We have found good agreement both for the shift and 
the damping with the experiment for m = 2, agreement 
also for the damping-rate of the m — mode, but dis- 
agreement similar to in the frequency shift for the 
m = mode. We explained this discrepancy with ob- 
servation by an artefact of our perturbative treatment 
and by applying the argument of Ref. |^,^ , namely by 
supposing an other nearby lying to = mode. 

To decide, whether a nonperturbative treatment of the 
problem, but still within the framework of our model ap- 
proximation, can show the existence of the missing m = 
mode requires an extended numerical study. Our future 
plan is to find the poles in a nonperturbative way and 
study the decoupling of the single particle and the col- 
lective excitations for T ^ Tc — 0, along the lines of 
the calculation for homogeneous system of Ref. ||ll| . A 
further natural extension of this work could be to study 
further and more difficult, but numerically still manage- 
able model approximations. 



ACKNOWLEDGMENTS 

This work has been supported by a project of the 
Hungarian Academy of Sciences and the Deutsche 
Forschungsgemeinschaft under Grant No. 95. R. G. 
and J. R. wish to acknowledge support by the Deutsche 
Forschungsgemeinschaft through the Sonderforschung- 
bereich 237 "Unordnung and grofie Fluktuationen" . Two 
of us (A. Cs. and P. Sz.) would like to acknowledge sup- 
port by the Hungarian National Scientific Research Foun- 
dation under Grant Nos. OTKA T017493, T025866 and 
F020094, and by the Ministry of Education of Hungary 
under Grant No. FKFP1059/1997. 



[1] E. A. Cornell, J. R. Ensher, C. E. Wiemann, 

cond-mat/9903109 . 
[2] W. Ketterle, D. S . Durfee and D. M. Stamper-Kurn 

cond-mat/9904034 . 
[3] D. S. Jin, J. R. Ensher, M. R. Matthews, C. E. Wiemann 

and E. A. Cornell, Phys. Rev. Lett. 77, 420 (1996). 
[4] D. S. Jin, M. R. Matthews, J. R. Ensher, C. E. Wiemann 

and E. A. Cornell, Phys. Rev. Lett. 78, 764 (1997). 
[5] M. O. Mewes, M. R. Andrews, N. J. van Drouten, D. M. 

Kurn, C. G. Townsend and W. Ketterle, Phys. Rev. Lett. 

77, 988 (1996). 
[6] D. M. Stamper-Kurn, H.-J. Miesner, S. Innouye, M. R. 

Andrews and W. Ketterle, Phys. Rev. Lett. 81, 500 

(1998). 

[7] A. Griffin, Excitations in a Bose-Condensed Liquid 
(Cambridge University Press, Cambridge 1993). 

[8] Shang-Keng Ma and Chia-Wei Woo, Phys. Rev. 159, 165 
(1967). 

[9] I. Kondor and P. Szepfalusy, Acta Phys. Acad. Sci. Hung. 
24, 81 (1968). 

[10] V. K. Wong and H. Gould, Ann. Phys. (N. Y.) 83, 252 
(1974). 

[11] P. Szepfalusy and I. Kondor, Ann. Phys. (N. Y.) 82, 1 
(1974). 

[12] S. H. Payne and A. Griffin, Phys. Rev. B32, 7199 (1985). 
[13] Gy Bene, P. Szepfalusy Phys. Rev. A58, R3391 (1998). 
[14] A. Minguzzi and M.P. Tosi, J. Phys; Cond. Matter 9, 
10211 (1997) 

[15] F. Dalvovo, S. Giorgini, L. Pitaevskii, and S. Stringari, 

Rev. Mod. Phys. 71, 463 (1999). 
[16] W. Vincent Liu, Phys. Rev. Lett. 79, 4056 (1997). 
[17] L. P. Pitaevskii and S. Stringari, Phys. Lett. A 235, 

(1997). 

[18] S. Giorgini, Phys. Rev. A57, 2949 (1998).. 

[19] P. O. Fedichev, G. V. Shlyapnikov and J. T. M. Wahaven, 

Phys. Rev. Lett. 80, 2269 (1998). 
[20] P. O. Fedichev and G. V. Shlyapnikov, Phys. Rev. A58, 

3146 (1998). 

[21] D. A. W. Hutchinson, R. J. Dodd and K. Burnett, Phys. 
Rev. Lett. 81, 2198 (1998). 



9 



[22] M. J. Bijlsma a nd H. T. C. Stoof, |coiid-mat/9902065 



[23] H. T. C. Stoof, cond-mat/981242 



[24] J. Reidl, A. Csordas, R. Graham, P. Szepfalusy, Phys. 

Rev. A59, 3816 (1999). 
[25] . . . , in preparation 

[26] M. R. Andrews, D. M. Kurn, H.-J. Miesner, D. S. Durfee, 
C. G. Townsend, S. Inouye and W. Ketterle, Phys. Rev. 
Lett. 79, 553 (1997) 

[27] M. Fhesser, private communication. 



FIG. 1. Spectral contributions to the thermally excited 
atoms ^{lo) as a function of the frequency uj. We com- 
pare the result of the local density approximation (open cir- 
cles) with the histograms from the direct diagonalization for 
T = 95 nK and Nc = 7500. 



FIG. 2. Damping rates of the m — (triangles) and m = 2 
(circles) modes as a function of T' (for the definition of T' see 
Ref. [Q). Experimental values (_full symbols) with error bars 
are taken from Fig. 3. of Ref. M. Open symbols denote the 
theoretical predictions ^uji of our model calculation using the 
dielectric formalism. 



FIG. 3. Excitation frequencies for the m = and m = 2 
modes as a function of T' . Notations for T' and for the ex- 
perimental values are the same as in Fig. ^. Open symbols 
denote our theoretical values + 5Ra;i in units of the axial 
trap frequency. 
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